---
#title: "Models Networks"
#author: "Laura Gómez Murillo"
#output: html_document

#All code ran in R version 4.1.2 (2021-11-01)

#uploading libraries and data

library(tidyverse)
library(dplyr)
library(plyr)
library(mgcv)
library(MuMIn)
library(betareg)
library(glmmTMB)
library(DHARMa)
library(cowplot)
library(stringr)

#####===========
#####===================Network metrics analyses ==================================#####
#####=============Creates Fig. 2 in paper

#analyze how rainfall and habitat suitability influence network size, mean normal degree, weighted degree, skewness, clustering

metrics <- read_csv('Networkmetrics_Gomez_Tarwater_Fig2.csv')

#scale the variables
s.metrics <- metrics %>% 
  dplyr::mutate(s.averain= scale(averain)) %>% 
  dplyr::mutate(s.perc_suitable= scale(perc_suitable)) %>% 
  dplyr::mutate(site=as.factor(site)) %>% 
  dplyr::mutate(month=as.factor(month)) %>% 
  dplyr::mutate(flock_id=as.factor(flock_id)) %>% 
  dplyr::mutate(swarm_id=as.factor(swarm_id))

#betas don't like 0s or 1s, so changing close to 0 and 1
proportions <- s.metrics %>% 
  dplyr::mutate(meanND=replace(meanND, which(meanND==0),0.0000000001)) %>% 
  dplyr::mutate(meanND=replace(meanND, which(meanND==1),0.9999999999)) %>%
  dplyr::mutate(cluster=replace(cluster, which(cluster==0),0.0000000001)) %>% 
  dplyr::mutate(cluster=replace(cluster, which(cluster==1),0.9999999999))


###==================== Network size =======================### 
#Testing different random effects first 
m1<-glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site),
                 data = s.metrics, family=poisson)
m2<-glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = s.metrics, family=poisson)
m3<- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|month),
                 data = s.metrics, family=poisson)
m4 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id),
                 data = s.metrics, family=poisson)
m5 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) + (1|month),
                 data = s.metrics, family=poisson)
m6 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|swarm_id) + (1|month),
                 data = s.metrics, family=poisson)
m7 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id) + (1|month),
                 data = s.metrics, family=poisson)
m8 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id),
                 data = s.metrics, family=poisson)
m9 <- glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id) + (1|month),
                 data = s.metrics, family=poisson)

summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)
summary(m7)
summary(m8)
summary(m9)
#swarmID is the only one that explains any variation, the rest explain 0. Model 2 is now the global model

# Testing if we should include quadratic effect or not
m.a<-glmmTMB(N_size ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = s.metrics, family=poisson)
m.b<-glmmTMB(N_size ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = s.metrics, family=poisson)

AIC(m.a, m.b) # model without quadratic effect is the final model 

#global model
m_size<-glmmTMB(N_size ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = s.metrics, family=poisson)
summary(m_size)
confint(m_size)

#model validation 
qqnorm(resid(m_size),pch=16); qqline(resid(m_size))
testDispersion(m_size)
simulationOutput <- simulateResiduals(fittedModel = m_size, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

E1 <- resid(m_size, type = "pearson")
N <- nrow(s.metrics)
p <- length(coef(m_size))
Dispersion <- sum(E1^2) / (N - p)
Dispersion

#everything looks good with the model, so now we can predict

#predictions for the model
s.averain=seq(from=min(s.metrics$s.averain), to=max(s.metrics$s.averain)*0.95,by=0.1)
s.perc_suitable= mean(s.metrics$s.perc_suitable)

newdata <- data.frame(expand.grid(swarm_id=unique(s.metrics$swarm_id), s.averain=s.averain,s.perc_suitable=s.perc_suitable))
newdata$fitlink <- predict(m_size, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(m_size, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- exp(newdata$fitlink)
newdata$upper.ci <- exp(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- exp(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.averain * attr(s.metrics$s.averain, 'scaled:scale') + attr(s.metrics$s.averain, 'scaled:center')
newdata 

#plot
NetSize <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "" , y = "Network size")+
  scale_x_continuous(limits=c(2100,2900), breaks= seq(2100, 2900, by = 100),expand=c(0,0)) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),axis.text.x=element_blank(),axis.text.y=element_text(size=14))
NetSize


###============== Mean normalized degree (proportion data)==================####

#Testing random effects first 
m1 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site),
                 data = proportions, family=beta_family)
m2 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = proportions, family=beta_family)
m3 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|month),
                 data = proportions, family=beta_family)
m4 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id),
               data = proportions, family=beta_family)
m5 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) + (1|month),
                 data = proportions, family=beta_family)
m6 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|swarm_id) + (1|month),
                 data = proportions, family=beta_family)
m7 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id) + (1|month),
               data = proportions, family=beta_family)
m8 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id),
                 data = proportions, family=beta_family)
m9 <- glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id) + (1|month),
                 data = proportions, family=beta_family)

summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)
summary(m7)
summary(m8)
summary(m9)
#swarmID is the main one that explains any variation, the rest have litte to 0 effect. Model 2 is now the global model

# Testing if we should include quadratic effect or not
m.a<-glmmTMB(meanND ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = proportions, family=beta_family)

m.b<-glmmTMB(meanND ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = proportions, family=beta_family)
AIC(m.a, m.b) # model without quadratic effect is better

# final model
m_meanND <-glmmTMB(meanND ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = proportions, family=beta_family)
summary(m_meanND)
confint(m_meanND) #nothing came out as important here, so don't need to graph

#model validation 
qqnorm(resid(m_meanND),pch=16); qqline(resid(m_meanND))
testDispersion(m_meanND)
simulationOutput <- simulateResiduals(fittedModel = m_meanND, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

E <- resid(m_meanND)
hist(E, xlab="Residuals", main="")
F3 <- fitted(m_meanND)
plot( F3, E)
plot(proportions$s.averain, E)
plot(proportions$s.perc_suitable, E)

###======================== Weigthed degree (continuos data)================####
# I checked before between Gaussian and Gamma distributions and based on DHARMa results, I choose Gaussian 

#checking random effects structure
m1 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site),
                 data = s.metrics , family=gaussian)
m2 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = s.metrics , family=gaussian)
m3 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|month),
                 data = s.metrics , family=gaussian)
m4 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id),
                 data = s.metrics , family=gaussian)
m5 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) + (1|month),
                 data = s.metrics , family=gaussian)
m6 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|swarm_id) + (1|month),
                 data = s.metrics , family=gaussian)
m7 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id) + (1|month),
                 data = s.metrics , family=gaussian)
m8 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id),
                 data = s.metrics , family=gaussian)
m9 <- glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id) + (1|month),
                 data = s.metrics , family=gaussian)


summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)
summary(m7)
summary(m8)
summary(m9)
#m2 with swarm id was the best

# Testing if we should include quadratic effect or not
m.a<-glmmTMB(meanWD ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = s.metrics , family=gaussian)
m.b<-glmmTMB(meanWD ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = s.metrics , family=gaussian)
AIC(m.a, m.b) # model without quadratic effect is better

# final model
m_meanWD <-glmmTMB(meanWD ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = s.metrics, family=gaussian)
summary(m_meanWD)
confint(m_meanWD)

#model validation
qqnorm(resid(m_meanWD),pch=16); qqline(resid(m_meanWD))
testDispersion(m_meanWD)
simulationOutput <- simulateResiduals(fittedModel = m_meanWD, plot = F)
residuals(simulationOutput)
plot(simulationOutput) #all looks good here

#### Predicting for rainfall
s.averain=seq(from=min(s.metrics$s.averain),to=max(s.metrics$s.averain),by=0.1)
s.perc_suitable= mean(s.metrics$s.perc_suitable)

newdata <- data.frame(expand.grid(s.averain=s.averain, s.perc_suitable=s.perc_suitable, swarm_id=unique(s.metrics$swarm_id)))
newdata$fitlink <- predict(m_meanWD, newdata,re.form=NA,se.fit = TRUE)$fit
newdata$selink <- predict(m_meanWD, newdata, re.form=NA,se.fit = TRUE)$se.fit
newdata$realval <- newdata$fitlink
newdata$upper.ci <- newdata$fitlink + 1.96*newdata$selink
newdata$lower.ci <- newdata$fitlink - 1.96*newdata$selink
newdata$origlength<-newdata$s.averain  * attr(s.metrics$s.averain , 'scaled:scale') + attr(s.metrics$s.averain , 'scaled:center')
newdata 

#plot
WDrain <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "", y = "Weighted degree")+ 
  scale_x_continuous(limits=c(2100,2900), breaks= seq(2100, 2900, by = 100),expand=c(0,0)) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),axis.text.x=element_blank(),axis.text.y=element_text(size=14))
WDrain

#now predicting for habitat suitability
s.averain=mean(s.metrics$s.averain)
s.perc_suitable=seq(from=min(s.metrics$s.perc_suitable),to=max(s.metrics$s.perc_suitable),by=0.1)

newdata <- data.frame(expand.grid(s.averain=s.averain, s.perc_suitable=s.perc_suitable, swarm_id=unique(s.metrics$swarm_id)))
newdata$fitlink <- predict(m_meanWD, newdata,re.form=NA,se.fit = TRUE)$fit
newdata$selink <- predict(m_meanWD, newdata, re.form=NA,se.fit = TRUE)$se.fit
newdata$realval <- newdata$fitlink
newdata$upper.ci <- newdata$fitlink + 1.96*newdata$selink
newdata$lower.ci <- newdata$fitlink - 1.96*newdata$selink
newdata$origlength<-newdata$s.perc_suitable  * attr(s.metrics$s.perc_suitable , 'scaled:scale') + attr(s.metrics$s.perc_suitable , 'scaled:center')
newdata 

#plot
WDhs <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "", y = "Weighted degree")+ 
  scale_x_continuous(limits=c(0.4,1.0), breaks= seq(0.4, 1.0, by = 0.1),expand=c(0,0))+ scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),axis.text.x=element_blank(),axis.text.y=element_text(size=14))
WDhs

###========================= Cluster (proportion data)=============================#

#Removing NAs from data set
proportions1 <- proportions %>% drop_na(cluster)
hist(proportions1$cluster)
boxplot(proportions1$cluster)

#Testing random effects first 
m1 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site),
                 data = proportions1, family=beta_family)
m2 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = proportions1, family=beta_family)
m3 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|month),
                 data = proportions1, family=beta_family)
m4 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id),
                 data = proportions1, family=beta_family)
m5 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) + (1|month),
                 data = proportions1, family=beta_family)
m6 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|swarm_id) + (1|month),
                 data = proportions1, family=beta_family)
m7 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id) + (1|month),
                 data = proportions1, family=beta_family)
m8 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id),
                 data = proportions1, family=beta_family)
m9 <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id) + (1|month),
                 data = proportions1, family=beta_family)
summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)
summary(m7)
summary(m8)
summary(m9)
#m2 with swarm id was the best

# checking quadratic effect
m.a <- glmmTMB(cluster ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = proportions1, family=beta_family)
m.b <- glmmTMB(cluster ~ s.averain +  s.perc_suitable +  (1|swarm_id),
                 data = proportions1, family=beta_family)
AIC(m.a, m.b) #without quadratic is the best


#final global model
m_cluster<-glmmTMB(cluster ~ s.averain +  s.perc_suitable +  (1|swarm_id),
                 data = proportions1, family=beta_family)
summary(m_cluster)
confint(m_cluster)

#model validation
E <- resid(m_cluster)
hist(E, xlab="Residuals", main="")
F3 <- fitted(m_cluster)
plot( F3, E)
plot(proportions1$s.averain, E)
plot(proportions1$s.perc_suitable, E)
qqnorm(resid(m_cluster),pch=16); qqline(resid(m_cluster))
testDispersion(m_cluster)
simulationOutput <- simulateResiduals(fittedModel = m_cluster, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#predicting for rainfall
s.averain=seq(from=min(proportions1$s.averain), to=max(proportions1$s.averain),by=0.1)
s.perc_suitable= mean(proportions1$s.perc_suitable)
newdata <- data.frame(expand.grid(swarm_id=unique(proportions1$swarm_id), month=unique(proportions1$month), site=unique(proportions1$site), s.averain=s.averain,s.perc_suitable=s.perc_suitable))
newdata$fitlink <- predict(m_cluster, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(m_cluster, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- plogis(newdata$fitlink)
newdata$upper.ci <- plogis(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- plogis(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.averain * attr(proportions1$s.averain, 'scaled:scale') + attr(proportions1$s.averain, 'scaled:center')
newdata 

#plot
NetCluster <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "Mean annual rainfall" , y = "Clustering")+
  scale_x_continuous(limits=c(2100,2900), breaks= seq(2100, 2900, by = 100),expand=c(0,0)) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
 theme(text = element_text(size = 14),axis.text.y=element_text(size=14),axis.text.x=element_text(size=14))

NetCluster
NetCluster <- NetCluster + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
NetCluster


###========================= Skew ND (Gaussian)=======================##
SkewData <- s.metrics %>% drop_na(skewND)
hist(SkewData$skewND)
boxplot(SkewData$skewND)

# I checked before between Gaussian and Gamma and choose Gaussian after DHARMa checking 
m1 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site), 
                 data = SkewData , family=gaussian)
m2 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = SkewData , family=gaussian)
m3 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|month),
                 data = SkewData , family=gaussian)
m4 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id),
                 data = SkewData , family=gaussian)
m5 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) + (1|month),
                 data = SkewData , family=gaussian)
m6 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|swarm_id) + (1|month),
                 data = SkewData , family=gaussian)
m7 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site) +  (1|swarm_id) + (1|month),
                 data = SkewData , family=gaussian)
m8 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id),
                 data = SkewData , family=gaussian)
m9 <- glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable + (1|site/swarm_id) + (1|month),
                 data = SkewData , family=gaussian)

summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)
summary(m7)
summary(m8)
summary(m9) #m2 with swarm id was the best

#checking quadratic effect
m.a<-glmmTMB(skewND ~ s.averain + I(s.averain^2)+ s.perc_suitable +  (1|swarm_id),
                 data = SkewData , family=gaussian)
m.b<-glmmTMB(skewND ~ s.averain +  s.perc_suitable +  (1|swarm_id),
                 data = SkewData , family=gaussian)
AIC(m.a, m.b) #without quadratic is better


#final global model
m_skew<-glmmTMB(skewND ~ s.averain + s.perc_suitable +  (1|swarm_id),
                 data = SkewData , family=gaussian)
summary(m_skew)
confint(m_skew)

#model validation
qqnorm(resid(m_skew),pch=16); qqline(resid(m_skew))
testDispersion(m_skew)
simulationOutput <- simulateResiduals(fittedModel = m_skew, plot = F)
residuals(simulationOutput)
plot(simulationOutput)#looks good

#predicting for rainfall
s.averain=seq(from=min(SkewData$s.averain),to=max(SkewData$s.averain),by=0.1)
s.perc_suitable= mean(SkewData$s.perc_suitable)

newdata <- data.frame(expand.grid(s.averain=s.averain, s.perc_suitable=s.perc_suitable, swarm_id=unique(SkewData$swarm_id)))
newdata$fitlink <- predict(m_skew, newdata,re.form=NA,se.fit = TRUE)$fit
newdata$selink <- predict(m_skew, newdata, re.form=NA,se.fit = TRUE)$se.fit
newdata$realval <- newdata$fitlink
newdata$upper.ci <- newdata$fitlink + 1.96*newdata$selink
newdata$lower.ci <- newdata$fitlink - 1.96*newdata$selink
newdata$origlength<-newdata$s.averain  * attr(SkewData$s.averain , 'scaled:scale') + attr(SkewData$s.averain , 'scaled:center')
newdata 

#plot
Skewrain <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "", y = "Skewness")+ 
  scale_x_continuous(limits=c(2100,2900), breaks= seq(2100, 2900, by = 100),expand=c(0,0)) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),axis.text.x=element_blank(),axis.text.y=element_text(size=14))
Skewrain

#now predicing for habitat suitability
s.averain=mean(SkewData$s.averain)
s.perc_suitable=seq(from=min(SkewData$s.perc_suitable),to=max(SkewData$s.perc_suitable),by=0.1)

newdata <- data.frame(expand.grid(s.averain=s.averain, s.perc_suitable=s.perc_suitable, swarm_id=unique(SkewData$swarm_id)))
newdata$fitlink <- predict(m_skew, newdata,re.form=NA,se.fit = TRUE)$fit
newdata$selink <- predict(m_skew, newdata, re.form=NA,se.fit = TRUE)$se.fit
newdata$realval <- newdata$fitlink
newdata$upper.ci <- newdata$fitlink + 1.96*newdata$selink
newdata$lower.ci <- newdata$fitlink - 1.96*newdata$selink
newdata$origlength<-newdata$s.perc_suitable  * attr(SkewData$s.perc_suitable , 'scaled:scale') + attr(SkewData$s.perc_suitable , 'scaled:center')
newdata 

#plot
Skewhs <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "Habitat suitability", y = "Skewness")+ 
  scale_x_continuous(limits=c(0.4,1.0), breaks= seq(0.4, 1.0, by = 0.1),expand=c(0,0))+ scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),axis.text.x=element_text(size=14),axis.text.y=element_text(size=14))
Skewhs
Skewhs <- Skewhs + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
Skewhs


#now graph them all together to create Fig 2
WDrain   <- WDrain   + theme(plot.margin = margin(t=5, r=5, b=5, l=2))
WDhs     <- WDhs     + theme(plot.margin = margin(t=5, r=5, b=5, l=2))
Skewrain <- Skewrain + theme(plot.margin = margin(t=5, r=5, b=5, l=2))
Skewhs   <- Skewhs   + theme(plot.margin = margin(t=5, r=5, b=5, l=2))
NetSize  <- NetSize  + theme(plot.margin = margin(t=5, r=5, b=5, l=2))
NetCluster <- NetCluster + theme(plot.margin = margin(t=5, r=5, b=5, l=2))


All <- plot_grid(
  WDrain, WDhs,
  Skewrain, Skewhs,
  NetSize, NULL,
  NetCluster, NULL,
  ncol = 2,
  nrow = 4,
  labels = c("a","b","c","d","e",NA,"f",NA),
  label_size = 12,
  label_fontface = "bold",
  label_x = 0.02,
  label_y = 1.02,
  hjust = 0,
  align = "hv",
  axis = "tblr",
  rel_widths = c(1,1),
  rel_heights = c(1,1,1,1)
)

save_plot("Fig2_Oecologia.png",
          All,
          base_width = 7,
          base_height = 9,
          dpi = 600)

#####===========
#####===================Dissimilarity analyses ==================================#####
#####=============
#analyze how rainfall and habitat suitability influence dissimilarity 

#### Create Fig 3 in the paper

#graphing the average of dissimilarity per site 
dissim.pairs <-read_csv('Networkoveralldissim_Gomez_Tarwater_Fig3.csv')

dissim.pairs$sites<-as.factor(dissim.pairs$sites)
dissim.pairs$dissimilarity<-as.factor(dissim.pairs$dissimilarity)
#change site names 
  site_key <- c(
  GATS = 2,
  LIR2 = 1,
  JUAN = 3,
  LIMB = 6,
  OGRD = 5,
  PLRD = 4, 
  SHER = 7
)

site_codes <- names(site_key)
dissim.pairs$site1 <- str_extract(
  dissim.pairs$sites,
  paste(site_codes, collapse = "|")
)
dissim.pairs$site2 <- str_extract(
  dissim.pairs$sites,
  paste0("(", paste(site_codes, collapse = "|"), ")$")
)
dissim.pairs$site1_num <- site_key[dissim.pairs$site1]
dissim.pairs$site2_num <- site_key[dissim.pairs$site2]

dissim.pairs$site_pair <-
  paste0("site", dissim.pairs$site1_num,
         "_site", dissim.pairs$site2_num)
         
dissim.pairs$site_low  <- pmin(dissim.pairs$site1_num,
                               dissim.pairs$site2_num)
dissim.pairs$site_high <- pmax(dissim.pairs$site1_num,
                               dissim.pairs$site2_num)

dissim.pairs$site_pair <-
  paste0("site", dissim.pairs$site_low,
         "_site", dissim.pairs$site_high)

dissim.pairs[, c("sites", "site1", "site2",
                 "site1_num", "site2_num", "site_pair")]
                 
#plot
dis.plot <- ggplot(dissim.pairs, aes(fill= dissimilarity, x= site_pair, y= value)) +
  geom_bar(position = "stack",stat= "identity") +
  scale_fill_manual(values= c("#9ebcda", "#8856a7"), labels=c('Rewiring (βos)', 'Turnover (βst)'))+
  theme_classic()+ 
  theme(axis.text.x =  element_text(size = 16, angle=90))+
  labs(y="Overall dissimilarity (βwn)", x= "")+ 
  scale_y_continuous(breaks= c(0.20, 0.40, 0.60, 0.80, 1.0))+
  theme(text = element_text(size = 16)) 
  
dis.plot
ggsave("Overalldissimilarity_Fig3.jpg", plot = dis.plot, width = 8, height = 6, dpi = 300)
  
#####==================Now look at rainfall and habitat suitability  
     
#upload data 
dissim1 <- read_csv('Networkdissim_Gomez_Tarwater_Fig4.csv')

dissim2 <- dissim1 %>% 
  filter(site1!=site2)

dissim3 <- dissim2 %>% 
  unite(sites, c("site1", "site2")) %>% 
  drop_na()
  
#scale variables  
s.dissim <- dissim3 %>% 
  dplyr::mutate(s.diffrain= scale(diffrain)) %>% 
  dplyr::mutate(s.diffhs= scale(diffhs)) %>% 
  dplyr::mutate(s.relST= scale(relST)) %>% 
  dplyr::mutate(s.relOS= scale(relOS))

s.dissim$repair<-as.factor(s.dissim$repair)

###========= WN - overall dissimilarity===================####
#model for overall
m_WN <- gamm(WN~ s(s.diffrain,k=3,bs="cs") + s(s.diffhs,k=3,bs="cs"), random=list(repair=~1), data = s.dissim, family=betar)
summary(m_WN$gam)
plot(m_WN$gam)

#predicting for rainfall
s.diffrain = seq(from =min(s.dissim$s.diffrain), to=max(s.dissim$s.diffrain), length = 100)
s.diffhs = mean(s.dissim$s.diffhs)

newdata <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)
newdata$fitlink <- predict(m_WN, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(m_WN, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- plogis(newdata$fitlink)
newdata$upper.ci <- plogis(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- plogis(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.diffrain * attr(s.dissim$s.diffrain, 'scaled:scale') + attr(s.dissim$s.diffrain, 'scaled:center')
newdata 

#plot
disWNrain <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "" , y = "Overall dissimilarity\n(βwn)")+
  scale_x_continuous(breaks = c(300,600,900,1200,1500),expand = expansion(mult = c(0, 0))) + scale_y_continuous( expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 12),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12)) 

# predicting for habitat suitability
s.diffrain = mean(s.dissim$s.diffrain)
s.diffhs = seq(from =min(s.dissim$s.diffhs), to=max(s.dissim$s.diffhs), length = 100)

newdata1 <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)
newdata1$fitlink <- predict(m_WN, newdata1,re.form=NA,se.fit = TRUE, type="link")$fit
newdata1$selink <- predict(m_WN, newdata1, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata1$realval <- plogis(newdata1$fitlink)
newdata1$upper.ci <- plogis(newdata1$fitlink + 1.96*newdata1$selink)
newdata1$lower.ci <- plogis(newdata1$fitlink - 1.96*newdata1$selink)
newdata1$origlength<- newdata1$s.diffhs * attr(s.dissim$s.diffhs, 'scaled:scale') + attr(s.dissim$s.diffhs, 'scaled:center')

disWNhs <- ggplot(data= newdata1, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "" , y = "")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(breaks= c(0.80, 0.82, 0.84, 0.86,0.88,0.90), expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 14),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12)) 

###========= ST - species turnover===================####
m_ST <- gamm(s.relST~ s(s.diffrain,k=3,bs="cs") + s(s.diffhs,k=3,bs="cs"), random=list(repair=~1), data = s.dissim, family=gaussian)
summary(m_ST$gam)

#predicting for rainfall
s.diffrain = seq(from =min(s.dissim$s.diffrain), to=max(s.dissim$s.diffrain), length = 100)
s.diffhs = mean(s.dissim$s.diffhs)

newdata2 <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)
newdata2$fitlink <- predict(m_ST, newdata2,re.form=NA,se.fit = TRUE)$fit
newdata2$selink <- predict(m_ST, newdata2, re.form=NA,se.fit = TRUE)$se.fit
newdata2$realval <-(newdata2$fitlink) * sd(s.dissim$relST) + mean(s.dissim$relST)
newdata2$upper.ci <- newdata2$realval + 1.96 * newdata2$selink * sd(s.dissim$relST)
newdata2$lower.ci <- newdata2$realval - 1.96 * newdata2$selink * sd(s.dissim$relST)
newdata2$origlength<-newdata2$s.diffrain * attr(s.dissim$s.diffrain, 'scaled:scale') + attr(s.dissim$s.diffrain, 'scaled:center')

#plot
disSTrain <- ggplot(data= newdata2, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "" , y = "Proportion of overall\ndissimilarity due to turnover\n(βst/βwn)")+
  scale_x_continuous(breaks = c(300,600,900,1200,1500), expand = expansion(mult = c(0,0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 12),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12))  

#predicting for habitat suitability
s.diffrain = mean(s.dissim$s.diffrain)
s.diffhs = seq(from =min(s.dissim$s.diffhs), to=max(s.dissim$s.diffhs), length = 100)

newdata3 <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)
newdata3$fitlink <- predict(m_ST, newdata3,re.form=NA,se.fit = TRUE)$fit
newdata3$selink <- predict(m_ST, newdata3, re.form=NA,se.fit = TRUE)$se.fit
newdata3$realval <-(newdata3$fitlink) * sd(s.dissim$relST) + mean(s.dissim$relST)
newdata3$upper.ci <- newdata3$realval + 1.96 * newdata3$selink * sd(s.dissim$relST)
newdata3$lower.ci <- newdata3$realval - 1.96 * newdata3$selink * sd(s.dissim$relST)
newdata3$origlength<- newdata3$s.diffhs * attr(s.dissim$s.diffhs, 'scaled:scale') + attr(s.dissim$s.diffhs, 'scaled:center')
 
#plot
disSThs <- ggplot(data= newdata3, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "" , y = "")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(breaks= c(0.60,0.65, 0.70, 0.75, 0.80, 0.85), expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 12),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12)) 



###========= OS - species rewiring===================####
m_OS <- gamm(s.relOS~ s(s.diffrain,k=3,bs="cs") + s(s.diffhs,k=3,bs="cs"), random=list(repair=~1), data = s.dissim, family=gaussian)
summary(m_OS$gam)

#predicting for rainfall
s.diffrain = seq(from =min(s.dissim$s.diffrain), to=max(s.dissim$s.diffrain), length = 100)
s.diffhs = mean(s.dissim$s.diffhs)

newdata4 <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)
newdata4$fitlink <- predict(m_OS, newdata4,re.form=NA,se.fit = TRUE, type="link")$fit
newdata4$selink <- predict(m_OS, newdata4, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata4$realval <-(newdata4$fitlink) * sd(s.dissim$relOS) + mean(s.dissim$relOS)
newdata4$upper.ci <- newdata4$realval + 1.96 * newdata4$selink * sd(s.dissim$relOS)
newdata4$lower.ci <- newdata4$realval - 1.96 * newdata4$selink * sd(s.dissim$relOS)
newdata4$origlength<-newdata4$s.diffrain * attr(s.dissim$s.diffrain, 'scaled:scale') + attr(s.dissim$s.diffrain, 'scaled:center')
newdata4 

#plot
disOSrain <- ggplot(data= newdata4, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "Difference in mean annual rainfall" , y = "Proportion of overall\ndissimilarity due to species rewiring\n(βos/βwn)")+
  scale_x_continuous(breaks = c(300,600,900,1200,1500), expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 12),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12)) 

#predicting for habitat suitability
s.diffrain = mean(s.dissim$s.diffrain)
s.diffhs = seq(from =min(s.dissim$s.diffhs), to=max(s.dissim$s.diffhs), length = 100)

newdata5 <- data.frame (s.diffrain=s.diffrain, s.diffhs=s.diffhs)

newdata5$fitlink <- predict(m_OS, newdata5,re.form=NA,se.fit = TRUE, type="link")$fit
newdata5$selink <- predict(m_OS, newdata5, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata5$realval <-(newdata5$fitlink) * sd(s.dissim$relOS) + mean(s.dissim$relOS)
newdata5$upper.ci <- newdata5$realval + 1.96 * newdata5$selink * sd(s.dissim$relOS)
newdata5$lower.ci <- newdata5$realval - 1.96 * newdata5$selink * sd(s.dissim$relOS)
newdata5$origlength<- newdata5$s.diffhs * attr(s.dissim$s.diffhs, 'scaled:scale') + attr(s.dissim$s.diffhs, 'scaled:center')
newdata5 

#plot
disOShs <- ggplot(data= newdata5, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "Difference in habitat suitability" , y = "")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(breaks= c(0.15, 0.20, 0.25, 0.30, 0.35, 0.40), expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#2D708EFF") +
  geom_line() +
  theme_classic()+
  theme(text = element_text(size = 12),
        axis.text.y=element_text(size=12),axis.text.x=element_text(size=12)) 

#now put them all together to make Fig 4
All <- plot_grid(disWNrain, disWNhs, disSTrain, disSThs, disOSrain,disOShs, ncol=2, nrow = 3, labels = c("a", "d", "b", "e", "c", "f"), label_size = 12, align = "v")
All

ggsave("All_Fig4.jpg", plot = All, width = 8, height = 10, dpi = 300)